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I.  INTRODUCTION 


Recent  advances  m  the  efficient  evaluation  of  integral  derivatives 
have  resulted  in  a  renewed  interest  in  the  analytical  evaluation  of  the  first 
and  second  derivatives  of  the  energy  with  respect  to  nuclear  coordinates  for 
Hartree-Fock  and  post-Hartree-Fock  wavefunctions.  Much  of  this  work  has  been 
pioneered  by  Pople  and  co-workers^  ^  who  developed  efficient  codes  both  for 
the  evaluation  of  the  second  derivatives  of  the  integrals  and  for  solving  the 
coupled-perturbed  Hartree-Fock  equations.  The  first  and  second  derivatives 
are  extremely  useful  tools  for  optimizing  structures,  characterizing 
stationary  points,  and  calculating  vibrational  frequencies,  as  well  as  for 
more  extensive  studies  of  the  Born-Oppenheimer  potential  energy  surface  away 
from  the  stationary  points.  The  last  several  years  have  also  seen  tremendous 
advances  made  in  CI°“*®  and  MCSCF^-^?  methodologies,  which  coupled  with  the 
interest  in  derivatives,  have  led  to  powerful  techniques  for  a  wide  range  of 
wavefunctions.  In  this  vein,  we  report  our  work  on  general,  mult ireference  Cl 
gradients  and  MCSCF  second  derivatives. 


Before  presenting  our  work,  it  is  appropriate  to  quickly  review  the 
previous  activity  in  the  area  of  Cl  gradients  and  MCSCF  second  derivatives. 

In  1980,  both  Krishnan,  et  al .  and  Brooks,  et  al .  presented  both  equations 
and  algorithms  for  the  analytical  computation  of  single-reference  UHF  or 
closed-shell  RHF  Cl  gradients.  Osamura,  et  al  .  extended  the  capabilities  to 
include  open-shell  RHF  reference  functions, ^  and  more  recently  published  the 
equations  for  general  multireference  Cl  gradients.  ®  Yamaguchi ,  et  al.^* 
reported  second  derivatives  for  two-configuration  SCF  (TCSCF)  wavefunctions 
and  Camp,  et  al .  ^  developed  the  methodology  and  implemented  second 
derivatives  for  CAS  MCSCF  wavefunctions.  Independently,  Pulay,  3  as  well  as 
Jorgensen  and  Simons, ^  have  provided  equations  for  MCSCF  second  derivatives. 
The  work  on  both  mult ireference  Cl  gradients  and  MCSCF  second  derivatives  is 
closely  related  since  a  major  step  in  both  types  of  calculation  is  the 
solution  of  the  coupled-perturbed  MCSCF  (CPMCSCF)  equations. 

In  this  paper  we  report  the  first  general  multireference  Cl  gradient 
cal culat ions .  We  have  extended  the  theory  developed  by  Osamura,  et  al.^®  to 
include  a  general  class  of  references  in  the  Cl.  This  extension  is  necessary 
for  the  calculation  of  gradients  for  a  commonly  employed  class  of  Cl 
wavefunctions  for  which  the  reference  configurations  are  selected  from  a 
generalized  CAS  MCSCF  wavefunct ion.  By  generalized  CAS  we  mean  that  there  is 
at  least  one  partially  occupied  orbital  subspace  where  the  energy  is  invariant 
to  rotations  of  the  orbitals  in  space.  The  simplest  example  of  this  type  of 
wavefunction  is  a  full  Cl  in  a  selected  subspace  of  orbitals.  The  length  of 
the  Cl  expansion  in  such  a  calculation  grows  rapidly  with  the  size  of  the 
active  space.  To  avoid  unreasonably  large  MCSCF  expansions,  it  is  desirable 
to  partition  the  space  of  chemically  active  orbitals  and  perform  CAS 
calculations  in  some  or  all  of  the  subspaces,  while  a  GVB  type  wavefunction 
might  be  used  in  the  remaining  subspaces.  A  MCSCF  calculation  where  all 
single  and  double  excitations  are  generated  from  one  subspace  into  another 
subspace  would  also  be  considered  a  generalized  CAS  wavefunction.  In  order  to 
calculate  the  gradient  of  a  selected  reference  Cl  wavefunction  in  which  the 
orbitals  are  obtained  from  a  generalized  CAS  wavefunction,  orbital  derivatives 
not  appearing  in  the  CPMCSCF  equations  are  needed  and  we  develop  the  machinery 
necessary  to  obtain  these  quantities. 


5 


We  also  present  the  equations  employed  in  our  MCSCF  second  derivative 
calculations.  Our  equations  are  an  extension  of  the  open-shell  SCF  second 
derivative  formulas  of  Osamura,  et  al .  2  These  equations  are  different  from 
those  presentedy  by  Pulay^  and  Jorgensen  and  Simons^  principally  in  the 
manner  in  which  the  derivative  overlap  contributions  are  handled.  We  expand 
upon  these  differences  later  in  the  text.  In  addition,  the  open-shell  second 
derivative  formulas  of  Osamura,  et  al.  2  are  reformulated  so  that  we  need  not 
contract  Coulomb  or  exchange  operators  after  the  CPMCSCF  equations  are  solved. 

The  CPMCSCF  equations  are  also  presented  in  a  more  compact  and  convenient 
manner  than  Osamura' s.  The  gradient-like  terms  appearing  in  the  CPMSCF  are 
defined  in  terras  of  modified  integrals  involving  derivative  overlap  terms,  as 
well  as  integrals  involving  derivative  atomic  orbitals.  This  formulation  is 
particularly  convenient  when  MCSCF  second  derivatives  are  desired. 

Finally,  we  present  the  results  of  sample  calculations  on  the  reaction, 

Be  +  H2  ♦  BeH2.  Here,  second  derivatives  obtained  at  the  MCSCF  level  are  used 
in  the  optimization  of  the  geometry  of  the  products  and  for  locating  the 
transition  state  at  the  multireference  Cl  level.  We  compare  the  results 
obtained  at  the  second-order  Cl  level  with  the  results  obtained  in  a  selected 
reference  Cl  (SRCI)  calculation. 


II.  COMPUTATIONAL  DETAILS 


A.  First  Derivatives 

The  energy  of  a  general  Cl  wavefunction 
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where  C  is  the  Cl  vector,  X  is  a  configuration  state  function  (CSF),  h^j  and 

gijki  are  one  and  two  electron  MO  integrals,  and  are  spin  coupling 

constants,  and  Hp^  is  an  element  of  the  Hamiltonian  matrix  in  our  CSF  basis. 

D- •  and  D,-  are  one-particle  and  two-particle  density  matrix  elements, 
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respectively.  The  first  derivative  of  the  Cl  energy,  with  respect  to  nuclear 
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coordinates,  is 
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In  this  equation  haj  and  are  derivatives  of  the  one-  and  two-electron 


atomic  orbital  integrals  transformed  to  the  molecular  orbital  basis.  The 
derivatives  of  the  molecular  orbital  expansion  coefficients  corresponding  to 
orbital  i,  are  expressed  in  terms  of  undifferentiated  coefficient  vectors 
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Lr^  is  a  Lagrangian  multiplier  defined  as 
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For  a  MCSCF  wavefunction 
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These  relations  are  simply  a  reflection  of  the  variational  conditions  on  the 
orbitals.  Not  all  of  the  elements  of  Uaj  are  independent.  By  differentiating 
the  orthonormality  condition 
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one  obtains 
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is  the  derivative  of  the  x,y  atomic  orbital  overlap  integral  with  respect 
to  nuclear  coordinate  a.  Uaj  can  then  be  expressed  as  a  sum  of  an 
antisymmetric  matrix  and  an  upper  triangular  matrix 
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and  T|j  is  defined  by 
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Thus,  the  derivative  of  the  MCSCF  energy  is 
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Before  proceeding  with  a  discussion  of  the  CPMCSCF  equations,  it  is 
useful  to  further  examine  the  contributions  to  the  first  derivative.  We  note 
that  it  is  possible  to  rewrite  the  first  derivative  for  a  Cl  wavefunction  as 
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This  equation  is  obtained  by  reordering  the  suras  in  Eq .  (3a)  as  follows: 
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Similarly,  the  first  derivative  of  the  MCSCF  energy  is 
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Partial-derivative  integrals  similar 
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the  MCSCF  second  derivatives. 


B.  Coupled  Perturbed  MCSCF  Equations 

The  derivative  of  the  molecular  orbitals  ufj >  which  are  needed  to  compute 

the  Cl  gradient,  are  obtained  by  solving  the  CPMCSCF  equations.  These 
equations  are  generated  by  requiring  that  the  wavefunction  satisfy  the  MCSCF 
variational  conditions  to  first  order  with  a  change  in  nuclear  geometry. 

Thus,  these  equations  result  from  requiring  that  the  derivatives  of  the 
orbital  and  Cl  stationary  conditions,  with  respect  to  a  nuclear  coordinate, 
vanish . 
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The  prime  sign  indicates  that  the  sum  only  runs  over  the  unique  orbital 
rotations  which  change  the  energy.  Gathering  the  terms  which  are  known  on  the 
right-hand  side  of  the  equation,  we  obtain  the  CPMCSCF  equations 
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In  these  expressions  the  superscript  "a"  in  G|:  and  (Gp1)3  is  used  to  denote 
that  derivative  AO  integrals  are  employed  in  the  construction  of  these 
quantities.  The  left-hand  side  of  these  equations  contains  explicitly  the 
Hessian  matrix  appearing  in  second  order  MCSCF  theory.  These  equations  are  in 
fact  similar  in  structure  to  the  second  order  MCSCF  equations,  and  can  thus  be 
most  efficiently  solved  by  expressing  the  Cl  variations  in  the  CSF  basis,  as 
noted  by  Lengsfield  and  Liu.^’^  These  equations  are  also  similar  to  those 
obtained  by  Osamura,  et  al.^G  However,  we  have  reordered  the  sums  appearing 
on  the  right-hand  side  of  this  equation  in  order  that  the  quantities  needed  to 
compute  the  MCSCF  second  derivatives  are  readily  available. 

C.  MCSCF  Second  Derivatives 

The  derivative  with  respect  to  a  nuclear  displacement  can  be  expanded  as 
follows : 


d _ ,3_  all  all 

db  v3b  ft 


r  i  3 
3b  3U 


3  c  n  a 
r  — E  _L 
P  3b  3C 


(22a) 


10 


(22b) 


'3b 


all  all 

+  5  5 


U 


ri  3U 


ri 


+  l 


P  3C 


The  second  derivative  of  our  MCSCF  wavefunction  is  obtained  by  applying  this 
operator  to  our  expression  for  the  first  derivative.  It  simplifies  matters  if 
we  make  use  of  both  expressions  (Eqs.  (13)  and  (13))  for  the  first  derivative. 
Our  second  derivative  expression  is  obtained  by  operating  on  Eq .  (13)  with  the 
first  two  terms  in  Eq.  (22)  and  then  operating  on  Eq.  (13)  with  the  last  term 
in  Eq.  (22). 

We  obtain 
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Note  that  only  the  first  two  terms  are  unique  to  MCSCF  second  derivatives. 
The  remaining  terms  are  equivalent  to  the  open-shell  formulas  of  Osamura,  et 
al.^2  It  is  also  important  to  note  that  almost  all  of  the  quantities 
involving  MO  or  Cl  derivatives  are  generated  by  setting  up  or  solving  the 
CPMCSCF  equations.  In  particular,  the  vector 
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can  be  stored  when  +  Giy  is  constructed.  Moreover,  the  Lagrangian 

multipliers  Laj  and  L^;-  are  used  to  obtain  the  gradient  terms  G|j  and  G-j  in 

the  CPMCSCF  equations.  The  trace  of  the  density  matrix  with  the  second 
derivative  integrals  is  calculated  by  transforming  the  density  matrix  elements 
to  the  AO  basis.  The  most  time  consuming  step  in  this  transformation  only 
requires  mn^  multiplications  where  m  is  the  number  of  active  orbitals  and  n  is 
number  of  basis  functions.  The  remaining  contributions  can  be  obtained  from 
the  product  of  one-particle  density  matrices. 

tb 

This 


•  •  •  •  t  T*"1 

The  only  term  which  requires  further  consideration  is  T'j 


term  formally  requires  one  to  contract  Coulomb  and  exchange  operators  (or  more 


efficiently  to  contract  Osamura’s  Y 


.Uc 


mm  j 


matrix)  with  the  solution  to  the 


CPMCSCF  equations,  U  ,  to  form  L^.j. 
noting  that 


However,  we  can  eliminate  this  step  by 


•I1  °£C  Lu.  T*.  -  “P  °f  L  .  T*.  Ub. 
f  T  ri  rx  r  i  ri  n  ri 


Y-  L..  (Ub  T3.  -  TS  Ub.)  . 

ij  ji  P  rp  pi  rp  pi 


(24) 


The  right-hand  side  of  this  equation  is  particularly  convenient  as  all  of 
these  terms  are  generated  when  we  set  up  the  CPMCSCF  equations.  Using  this 
relation  and  Eq.  (10),  we  now  obtain  the  final  second  derivative  expression. 


rab 


a  T3 
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pq  pq 
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Equation  (23)  can  be  decomposed  in  another  fashion  if  the  orthonormality 
conditions  are  expressed  in  such  a  way  that  T  is  symmetric.  Camp,  King, 
Mclver,  and  Mullally  ^  have  expressed  the  orthonormality  conditions  in  this 
way.  However,  there  is  an  advantage  in  employing  an  upper  triangular  T  matrix 
as  only  the  first  p  (p  is  the  number  of  occupied  orbitals)  columns  are  needed 

to  obtain  the  integrals  hi.  and  Thus,  the  transformations  needed  to 

obtain  h-j  and  gjj^i  can  be  performed  very  efficiently.  This  fact  is  also 

exploited  in  Osamura's  equations. 

D.  Selected  Reference  Cl  Gradients 


The  gradient  of  the  energy  for  a  Cl  wavefunction  requires  knowledge  of 
the  first  order  variations  in  the  molecular  orbitals.  The  CPMCSCF  scheme  does 
not  uniquely  define  a  tranformat ion  of  the  orbitals,  but  only  specifies  a 
transformation  of  the  variational  parameters.  This  is  enough  information  to 
define  an  orbital  transformation  excluding  an  arbitrary  orthogonal  mixing 
within  the  invariant  subspaces.  If  the  Cl  wavefunction  has  the  same  invariant 
subspaces  as  the  reference  wavefunction,  then  the  gradient  is  well  defined. 
This  is  the  case,  e.g.,  with  a  SDCI  using  closed-shell  Hartree-Fock  orbitals. 
Both  wavefunctions  are  invariant  to  mixings  among  the  doubly  occupied  core 
orbitals  and  information  concerning  such  mixings  is  not  required  for  the  Cl 
gradient . 

For  a  Cl  wavefunction  constructed  as  all  excitations  of  a  given  order 
from  selected  references  of  a  generalized  CAS  wavefunction,  this  is  not  the 
case.  The  Cl  is  not  invariant  to  mixings  of  orbitals  within  the  partially 
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occupied  subspaces.  The  most  common  way  to  address  this  problem  in 
calculating  the  Cl  energy  itself  is  to  require  that  the  orbitals  be  natural 
orbitals  of  the  one-particle  density  matrix.  For  the  calculation  of  the  Cl 
gradient,  we  then  remove  the  ambiguity  by  requiring  the  natural  orbital 
conditions  to  be  satisfied  to  first  order  with  a  change  in  nuclear  geometry. 
The  derivative  of  the  natural  orbital  condition  must  vanish  in  analogy  with 
Eqs.  (18)  and  (19)  for  the  variational  conditions. 


■n  <D“  - w)  ■ 


(26) 


Here  D  is  a  subblock  of  the  one  particle  density  matrix  and  \  is  a 
diagonal  matrix  of  occupation  numbers.  At  a*“0,  W*I,  and  thus  D(0)=X(0). 


Evaluating  Eq.  (26)  leads  to  expressions  for  the  derivatives  of 
parameters  which  mix  orbitals  within  the  invariant  subspace. 


a  ,.a 
A  •  •  =  W  .  . 
iJ  ij 


-D3. 
_ LL_ 

Xi‘Xj 


(27) 


In  this  equation,  D|j  =  2  I  CpKpqDq  an<*  Cp  *s  obtained  from  the 
solution  to  the  CPMCSCF  equations.  Afj  is  then  obtained  with  the  overlap 
derivative  portion  T|j  as  in  Eq.  (10)  to  obtain  Ua ^ . 


III.  COMPARISON  WITH  OTHER  WORK 

As  noted  in  the  Introduction,  our  second  derivative  expressions  differ 
from  those  of  Pulay^  and  Jorgensen  and  Simons^  principally  in  the  manner  in 
which  the  derivative  overlap  terras  are  included  in  the  CPMCSCF  equations,  and 
subsequently  in  the  final  expression  for  the  second  derivatives.  Jorgensen 
and  Simons  include  an  overlap  term  in  the  definition  of  their  AO  basis.  They 
are  then  able  to  derive  a  very  compact  set  of  formulas.  However,  the  formulas 
for  the  derivatives  of  their  AO  integrals  are  involved  and  the  overall 
efficiency  of  their  method  depends  very  strongly  on  how  their  derivative  AOs 
are  computed.  We  should  also  note  that  Jorgensen  and  Simons'  method  was 
developed  to  describe  second  derivatives  of  Cl  and  coupled-cluster 
wavefunctions,  as  well  as  MCSCF  wavefunct ions .  The  comparative  efficiency  of 
their  method  must  also  be  judged  on  how  well  it  treats  highly  correlated  wave- 
functions,  and  this  analysis  is  beyond  the  scope  of  the  present  study. 

Pulay  has  derived  the  CPMCSCF  equations  in  such  a  manner  that  the  orbital 
variations  are  expressed  in  the  AO  basis  (as  opposed  to  the  MO  basis  used  by 
Osamura)  and  the  Cl  variations  are  in  the  CSF  basis.  His  CPMCSCF  equations  do 
not  require  a  contraction  of  Coulomb  and  and  exchange  operators  with  overlap 
terms,  but  his  CPMCSCF  equations  are  larger  than  Osamura' s  because  he  is 
working  in  the  AO  basis.  We  feel  that  the  overall  efficiency  of  both  methods 
should  be  about  the  same.  However,  we  also  note  that  Pulay' s  equations 
neglect  projection  operator  terms  appearing  in  Osaraura's  equations  (see  also 
Lengsfield  and  Liu^).  These  terras  are  needed  if  the  Hessian  appearing  in  the 
CPMCSCF  equations  is  to  be  nonsingular. 


The  recent  communication  by  Camp,  et  al . ,  *  provides  few  details  of  their 
second  derivative  equations.  However,  they  express  the  variation  of  the 
molecular  orbitals  with  nuclear  displacement  in  terms  of  a  product  of  an 
exponential  and  a  Hermitian  operator.  The  orthonormality  conditions  at  first 
order  are  treated  by  this  Hermitian  operator  as  opposed  to  the  upper 
triangular  matrix  appearing  in  Eq.  (12).  Further  details  are  needed  to 
determine  if  the  simplicity  of  their  final  equations  offset  the  computational 
expense  of  working  with  a  Hermitian  matrix  as  opposed  to  an  upper  triangular 

Ta 

matrix  in  constructing  terms  like  L  appearing  in  our  second  derivative 
expression. 


IV.  SAMPLE  CALCULATIONS 


The  reaction  Be  +  H^  *  BeH2  was  studied  in  Cjv  symmetry.  The  geometry  of 
the  products  and  transition  state  was  stabilized  at  the  MCSCF  and  multi¬ 
reference  Cl  level.  The  MCSCF  wavefunction  employed  in  this  study  was  a  four- 
electron  in  four-orbital  CAS.  This  wavefunction  correctly  describes  the 
cleavage  of  the  BeH  bonds,  but  does  not  contain  all  of  the  configurations 
needed  to  describe  s^  and  p^  near  degeneracy  in  Be.  Second  order  Cl 
calculations,  based  on  the  4  in  A  CAS,  are  compared  to  the  results  of  selected 
reference  singles  and  doubles  Cl  (SRCI)  calculations.  In  the  SRCI 
calculations,  the  references  were  selected  on  the  basis  of  their  cumulative 
weight  in  the  natural  orbital  representation  of  the  MCSCF  wavefunction.  The 
weight,  W,  was  defined  as  follows: 


W 


i  C! 


(28) 


OQ 

The  basis  set  used  for  beryllium  was  Dunning's  5s  contraction  °  of  Huzinaga's 
10s  primitive  set^  and  Bartlett's  p  function^®  composed  of  three  primitives 
which  we  augmented  with  an  uncontracted  p  function  of  exponent  0.057  181.  Th^ 
hydrogen  basis  was  Dunning's  2s  contraction^  of  Huzinaga's  4s  primitive  set"*' 
with  a  scale  factor  of  1.2.  The  results  of  these  calculations  are  presented 
in  the  following  four  tables.  Table  1  lists  the  references  employed  in  the 
SRCI  calculations.  The  stable  geometries  are  given  in  Table  2,  and  the  MCSCF 
vibrational  frequencies  in  Table  3.  Our  calculated  MCSCF  and  Cl  energies  are 
given  in  Table  4. 


The  reaction  is  symmetry  forbidden  along  the  C2V  reaction  path  and  this 
results  in  a  transition  state  with  two  dominant  CSFs. 


The  results  of  the  SRCI  calculations  accurately  reproduce  the  second 
order  Cl  results  while  only  requiring  a  fraction  of  the  computational  effort. 
The  MCSCF  calculations  provide  a  good  description  of  the  reaction  (in  this 
basis)  but  the  activation  energy  is  a  bit  high  as  expected. 

We  found  that  the  multireference  Cl  gradient  calculations  converged  very 
rapidly  to  stable  points  when  MCSCF  second  derivatives  and  starting  geometries 
were  employed.  In  the  most  unfavorable  case,  four  iterations  were  required 
for  convergence  (largest  component  of  the  gradient  less  than  1.0  x  10”  a.u ./ 
bohr) . 


Table  1.  Dominant  Configurations  in  the  MCSCF  Wavefunction 


Equilibrium  Geometry  Configuration 

1  la^ajlb* 

2  1 a^2a ^ 1 b2 3a ^ 2b2 

Weight  *  0.98766 

Transition  State  Geometry  Configuration 

1  la^a^lb^ 

2  Iaj2af3af 

3  Ia^2ajlb23aj2b2* 


Coefficient 

0.989 

0.092 

Coefficient 

0.740 

-0.623 

-0.197 

0.117 


♦The  two  coefficients  associated  with  this  configuration  correspond  to  the 
two  spin  couplings. 


Equilibrium 

RBeH 

0 

C2 

z\) 

0 


Table  2.  Geometries  of  BeH2  Stationary  Points 


MCSCF 


SRCI 


SDCI 


2.61  876  bohr 
180.00° 

3.0875  bohr 
47.38° 


2.60  512  bohr 
180.00° 

3.0687  bohr 
48.68' 


2.60  637  bohr 
180.00' 

3.0693  bohr 
48.76’ 


Table  3.  MCSCF  Vibrational  Frequencies 


Equilibrium 

790.1  cm"1 
790.1  cm"1 
1827.5  cm"1 
2069.4  cm"1 


Transition  State 
4044. 9 i  cm  1 

993.4  cm"1 
4514.7  cm"1 
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Table  4.  MCSCF  and  Cl  Energies  for  BeH2 

MCSCF  SRCI  SDCI 

Equilibrium  -15.773  651  a.u.  -15.798  597  a.u.  -15.798  965  a.u. 

Transition  State  -15.597  013  a.u.  -15.628  925  a.u.  -15.629  315  a.u. 

AE  0.176  638  a.u.  0.169  672  a.u.  0.169  650  a.u. 

(110.8  kcal)  (106.5  kcal)  (106.45  kcal) 


V.  CONCLUSION 

We  have  presented  a  set  of  simple  and  efficient  formulas  for  the 
calculation  of  multireference  Cl  gradients  and  MCSCF  second  derivatives.  The 
additional  machinery  needed  to  compute  selected-reference  Cl  gradients  was 
developed.  Sample  calculations  were  presented  in  which  SRCI  structures  and 
activation  energies  compared  very  favorably  with  the  full  second  order  Cl 
results.  In  the  reaction  investigated  in  this  work,  MCSCF  second  derivatives 
were  found  to  be  very  useful  in  the  stabilization  of  the  Cl  structures. 
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